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The path integral formalism gives a very illustrative and 
intuitive understanding of quantum mechanics but due to its 
difficult sum over phases one usually prefers Schrodinger's 
approach. We will show that it is possible to calculate simple 
quantum phenomena by performing Feynman's sum over all 
paths staying entirely in real time. Once the propagator is 
obtained it is particularly easy to get the energy spectrum or 
the evolution of any wavefunction. 

I. INTRODUCTION 

In 1942 R.P. Feynman presented the path integral 
formalism in configuration space as an alternative and 
very illustrative description of quantum mechanics [1,2]. 
Many applications in areas like quantum statistical me- 
chanics or field theory were explored in the following 
years [3-7]. 

Feynman described the probabilistic character of quan- 
tum mechanics by the ability of a quantum particle to 
take virtually all possible paths between two points. The 
different paths were weighted by a phase factor, depend- 
ing on the classical action of the path measured in units of 
Ti. The sum over all paths (the path integration) of these 
phase factors results in the amplitude for the probability 
of propagation. In the classical limit of very large ac- 
tions (in units of Ti) the Hamilton principle of stationary 
action becomes apparent, as all other highly oscillating 
phase factors can interfere destructively. 

But if we treat the genuine quantum mechanical case 
there arise severe mathematical difficulties due to the 
oscillating path phases. Many sophisticated methods 
have been developed to find the most important class 
of paths, e.g. stationary phase approximations (in semi- 
classical quantum mechanics) [13,8,9], or to damp the 
integrand (e.g. Euclidean (imaginary) time path integra- 
tion) as used for finite temperature calculations [3,4,6]. 

In the present work we will explore the direct use of 
the basic concept of the sum over all paths in configu- 
ration space and real time. We will perform it without 
any weighting or damping change of the phase factors in 
our numerical calculations by iterating the 1-time-slice 
propagator. It will be seen that simple quantum systems 
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such as the harmonic oscillator or the double well poten- 
tial can be well described by this illustrative method. 

The fact that one can perform the path integral di- 
rectly in Minkowski space could turn out to be important 
for many applications like for real time tunnelling effects 
at finite temperature. 

II. CALCULATION OF THE REAL-TIME 
PROPAGATOR 

A. Discretization method 

If we want to investigate the dynamics of a quantum 
system in configuration space we have to derive first the 
local representation of the time evolution operator U . 
Then, by using 

*(a;,tb) = < «|*(ib) > 

= j Ax' <x\U{ty„t^)\x' >-^{x' ,t^) (1) 

we obtain for any given initial wavefunction "^{x,t^) the 
propagated one "^{x,t\j) at a future time tb > ^a- 
The Feynman path integral, 

< «b| U (^b,^a)ka > 

f ^ t f*^ dtC(x,x) 

= Vx , (2) 

enables us to get this propagator, < ^bjt^l^^a >, from the 
sum over all paths between (a;a,^a) and (a;b,^b), symbol- 
ised by Vx. The C denotes the classical Lagrangian of 
the system. 

The fundamental composition law for the evolution op- 
erator suggests a discrete representation such as 

< «b| U (tb,^a)ka > 

= j Ax{ti)..Ax{tiq-i) 

N 

X ]^ < x{t„)\U {t„,t„_i)\x{t„_i) > 

n = l 

~ lim I dx{ti)...dx{tN-i) 

X exp 1^ '^^S(x(t„),t„,x(t„-i),t„-i)^ , (3) 

with the equidistant discretization (ti),i = 0..N, of the 
time interval [ta,^b] and (t^ = to, ty^ = t^^ , x^ = x(to), 
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xi, = x{tN), T := tb-ta, At := ^ and C := (^^^ffxt)")- 
The evaluation of the action value S needs care in non 
trivial geometries [3-5] but also in numerical calculations 
with finite N . We use the midpomi rule which can be seen 
in Eq. 4 below. 

Furthermore, we approximate the integrations over 
space by a finite equidistant discretization (xi), 
i = 0..-D, on a spatial interval somewhat larger than the 
region where the wave functions could be of importance, 
[xq, xo]- We will treat only time independent potentials 
which causes all propagators to be functions of just time 
differences in our examples. In order to get the prop- 
agator we break it up into short time intervals. Thus 
we can use as a good approximation for the short time 
propagator the matrix IC(At) defined by 

/C,,(AO := Cexp | — £ (^^, ^ j | 

« < Xi{U + At)\U{At)\xj{U) >, (4) 

for arbitrary values of Xi and Xj . We obtain then a uni- 
versal {D -\- l)-dimensional, complex matrix /C(At) inde- 
pendent of tk . 

Consequently, the approximated finite time propagator 

GijiT) « < Xi\U{T)\xj > 

results by taking the product over short time propagators 
IC{At) as 

GiT) := (Axf-' IC^iAt) , (5) 

with Ax := ^r>~^" ^ via simple matrix multiplications. 
The numerical effort can be substantially reduced by 
recursive matrix multiplication, needing therefore only 
ln(#/2) steps. Thus we can get now fairly precisely (with 
very small At) the probability amplitude of a quantum 
particle to go from Xi to Xj in time T . 

Interesting to compute is the propagation of a dis- 
cretized wavefunction ^ by inserting the approxima- 
tion (5) into (1). With it we can measure the evolving 
expectation values, 

<0>{t) :=< ^U+{t) 6 U(t)\^! > 

« J2 Shii) < ^j\0\^k > Gki(t) "^(xi), (6) 

for an arbitrary observable O in order to get physical 
information about its time evolution (e.g. the particle 
position or the conservation of energy). 

B. Energy Levels 



complete set of energy eigenstates \n>, with n = 0, 1, ... 
and 

H\n>=En\n>, (7) 
we can decompose 

Tr U{t) = j dx < x\e-^"*\x > 

= dx < x\e~i^*\n >< n\x > 

= ^e-i^-.* (8) 

n 

due to orthonormality of the wavefunctions. We obtain 
the energy levels after Fourier transformation^ of the 
trace [11]. However, even in analytic calculations this 
yields serious trouble because of singularities and poles, 
and needs to be addressed with care. 

The time development of the trace on the interval 
[T, (Nt + 1)T] can be performed recursively with the com- 
position law, 

G(ti) = G(T) G(ti-i) Ax, (9) 
for the propagator matrix (i = I..Nt and to = T). 

C. Normalization 

An estimate for the quality of the propagator can 
be obtained from the normalization condition, meaning, 
that the probability to go from one point to any possible 
point at a later time has to be 1, due to 

= /dA,.",<.T'i.'>r<«"ifi.>^ ,10, 

In our discrete formulation we should get 
D 

J2i^xf 6W)6^kit) = l (11) 

i,j=0 

for every starting point x^. 

In our examples we will see that this normalization is 
quite good in the needed spatial intervals. We will only 
use it to renormalize the long time propagator needed for 
the energy spectrum (Eq. 9) in order to compensate for 
numerical errors. 



One of the most important quantities to extract, for 
time independent systems, is the energy spectrum of the ^ We use the Fast Founer Transformation [12]. 
Hamiltonian H . We can get this from the propagator 
by calculating its trace as a function of time. Taking a 
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III. HARMONIC OSCILLATOR 

The harmonic oscillator represents the simplest exam- 
ple where it is possible to confine the calculation to a 
finite position interval. Furthermore numerical results 
can be compared with analytic solutions. 

In the numerical calculation we will use the scaled, one 
dimensional, classical Lagrangian C, 

C{x, x) = —x^ — 2*^- (12) 

with the resulting eigenfrequency To = 27r. For a suf- 
ficient time resolution of an oscillation we will perform 
the propagator Qij in intervals T = with 4 time slices 
(i.e. N = A and At = |^ in Eq. 5, respectively). 

First we investigate the time evolution for a real Gaus- 
sian wavefunction "^^{x), 

*o(a;) = (^)' e- = ("-"^*"*)', (13) 

which should not change as a generalized coherent state 
its starting shape after a period. It can be regarded as a 
quantum particle aX < x > which oscillates with a period 
Tq. This will be confirmed up to an error of < 10~® % if 
we measure the position expectation values successively 
between the propagations (Fig. 1). 



depends on the shape of the wavefunction. This can be 
seen in Fig. 2, where the propagation of the probabil- 
ity density |^(a;,nT)p for the above example (top) is 
shown. It is compared with a wavefunction which has 
a more classical-like time evolution of the energies (bot- 
tom). All numerical results are practically equal to the 
analytic ones. 
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FIG. 1. Energy and position expectation values as a func- 
tion of time for the harmonic oscillator (a = 2.0). 

Also is shown from measurements of the Hamiltonian 
that the conservation of energy is ensured up to < 10~® % 
per period. In addition, the virial theorem (e.g. in [10]) 
is well reproduced by estimating the time average of the 
potential energy. 

On looking more carefully at Fig. 1 a difference be- 
tween a classical and quantum particle can be extracted 
because the potential energy in quantum mechanics not 
necessarily possesses a maximum at the classical turn- 
ing points. The ratio of the kinetic and potential energy 
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FIG. 2. Time development of the probability density of a 
Gaussian wavefunction in the harmonic oscillator, with 
a = 2.0 (top) and a = 0.5 (bottom). 

For the calculation of the energy spectrum from the 
propagator trace (Eq. 8) we constructed a kernel with 
T = J. The energy range of the spectrum is determined 
by this value only. It is not affected by estimating the 
trace evolution for a longer time period but this increases 
the accuracy of the spectrum. The peaks appear better 
localized and give more precise energy eigenvalues. The 
resulting graph from a discrete Fourier transformation of 
the trace in time is shown in Fig. 3. Here the evolu- 
tion is calculated for a period of more than 70 harmonic 
oscillations. 
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The energy levels agree very well with the well known 
eigenvalues E„, 

Er, = n+^, n = 0,l,..., (14) 

of the harmonic oscillator Hamiltonian. So in this case 
the main physical characteristics can be reproduced with 
high precision by our realization of the Feynman path 
integral directly in real time. 



lution (e.g. in [2,3,11]) on the same lattice we get an 
estimate for the rather high quality of the norm. It elu- 
cidates that the use of a finite spatial interval causes the 
main deviations from the perfect norm. Notice that we 
can always choose the interval wider than the relevant 
allowed region for a given potential thus avoiding this 
problem. In addition the norm provides a helpful con- 
vergence criterion for the choice of the spatial grid. In 
our examples 600 grid points are sufficient to minimize 
the infiuence of the discretization in position space. 



IV. DOUBLE WELL POTENTIAL 
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FIG. 3. Energy spectrum of the harmonic oscillator from 
Fourier transformation of the propagator trace (Nt = 1023). 

Finally we show the propagator norm as a function of 
the starting position in Eq. 11 in Fig. 4. 
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FIG. 4. Norm of a harmonic oscillator propagator com- 
pared to the exact solution on the interval [—7, 7] for T = j^. 

Only after comparison with the exact continuous^ so- 



In a second example we examine the double well po- 
tential, 

jC{x, x) = - a{x - XMm f (x + XUm f , (15) 

having two harmonic oscillator like minima. This Mex- 
ican hat potential is important in quite different fields 
(cf. e.g. [3]) but there exists no closed analytic solution 
for the propagator. 

Here we want to show the time evolution of a tunnelling 
process as one of the most important phenomena in quan- 
tum mechanics. We construct an almost localized parti- 
cle in one well, from a superposition of the lowest lying 
symmetric and antisymmetric energy eigenstates and 
of the total system. They can be approximated with 
two parametric testfunctions Tg/A (^s '+' and Ta '— '); 
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/27r 



± e" 



(16) 



by a simple variational method. The displaced Gaussians 
in Eq. 16 correspond to a harmonic approximation in the 
two wells. The resulting superposition is shown in Fig. 5. 
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^In our example the discrete solution (cf. [3]) makes no sig- 
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FIG. 5. Construction of a localized wavefunction on one 
side from an approximation for the lowest lying symmetric 
(^'s) and antisymmetric (^'a) energy states of the double well 
potential. 
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Both energy levels are below the threshold in the po- 
tential centre and so a classical particle with the same en- 
ergy could not leave one of the wells. But this is possible 
in the quantum case as can be shown using our calculated 
propagator (Fig. 6). Here every 55th time step of the 
squared wavefunction is plotted (i.e. '■= U^^(T)^!o). 
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FIG. 6. Tunnelling of an initially (^'o) on one side localized 
wavefunction in a double well potential. 

Decisive for the tunnelling rate is the size of the barrier 
which determines the also drawn energy difference of the 
two constituents Ea — Eg (e.g. in the limit of very high 
barriers the energy levels become equal/degenerated). In 
our example the measurement of the position expectation 
value evolution (Fig. 7) results in a tunnelling time of 
r = 54.4. Furthermore the energy conservation is granted 
up to 1.2 % per tunnelling period. 
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7. Energy and position expectation value of a tun- 
wavefunction in a double well potential. 



The tunnelling time can be checked using 

U(t)( |*s > +|*A >) 

= e'^*(|*s>+|*A >) 



(17) 



assuming that we use the exact eigenstates. It is obvious 
that in a time tq := ^ the sign in the superpo- 

sition changes. This can be interpreted as a complete 
tunnelling to the other well. Precise measurements of 
the energy values Es and E^ in our numerical example 
gives a tunnelling time of tq = 61.4. This 13 % deviation 
should result partly from the inaccurate approximation 
of the energy states ^s/A- This can be proved by a look 
at the energy spectrum of the double well (Fig. 8). 
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FIG. 8. Energy spectrum of the double well potential from 
Fourier transformation of the trace. The potential barrier is 
low enough to show the lowest lying eigenvalues Es and i?A 
well separated(7VT 



2047, T = f ) 



Easy to see is the broken degeneracy of the lowest ly- 
ing eigenvalues Es and E^ due to the finite barrier. The 
peaks occur at 0.433 and 0.494, respectively, close to our 
Gaussian superposition {Es = 0.438, E^ = 0.490). How- 
ever, the energy difference in the eigenspectrum results 
in a tunnelling period of t^s = 51.7 which is only 5.0 % 
below our observed value r = 54.4 in Fig. 7. In addition 
the subsequent energy gaps are increasing as expected for 
all symmetric potentials which are of higher order than 
the harmonic one. 

Finally the normalization shows also the correct shape 
in accordance with our experience from the harmonic 
oscillator (Fig. 9). In our calculations coarser or non- 
equidistant spatial matrices make it worse. But on the 
other hand a finer grid hardly improves the quality and 
the computational effort remains reasonable. 
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FIG. 9. Norm of a propagator in a double well potential on 
the interval [-7,7] (T = f^). 

Both examples show that it is possible to describe sim- 
ple quantum phenomena with the illustrative Feynman 
real time path integral. It could be used not only as a 
didactic approach for the study of the elementary prop- 
erties of quantum theories in cases where one needs the 
real time evolution, which it is hard to get by analytical 
continuation. 
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